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ABSTRACT 


Future aircraft will eventually feature nonaxi symmetric or rectangular 
nozzles. Developing a three-dimensional code to simulate the characteristics 
of the jet exhaust plume, issuing from nonaxi symmetric nozzles, in general, at 
different flight conditions, is very important. In the present investigation, 
two three-dimensional codes were developed to simulate the shock-cell 
structure of circular nozzles. These codes are used to solve the parabolized 
and simplified Navier-Stokes equations respectively. Both codes are based on 
a method previously developed by Newsome et al. (Ref. 1). These codes are 
fully vectorized on the VPS 32 at NASA Langley Research Center. The 
axisymmetric underexpanded supersonic jet flow problem, exhausting into still 
air, was used as a test case for developing an efficient three-dimensional 
code which should be capable of simulating two-d'imensional problems and 
preserving crossplane symmetry of the flow downstream of the jet exit. 

INTRODUCTION 

Future propulsion systems for fighter aircraft must be designed for 
maximum maneuverability over a wide range of flight Mach numbers (Ref. 2 and 
3). These features can be achieved with rectangular or nonaxi symetric nozzles 
(Ref. 4-8). Rectangular nozzles provide more rapid plume velocity decay than 
axisymetric nozzles (circular jet) and is simpler to modify, for example, to 
incorporate thrust vectoring (Ref. 9). Developing an efficient computational 
technique is essential to fully understanding the flow characteristics of 
these nozzles (rectangular nozzles). This computational technique should have 
a three-dimensional calculation capability and be able to simulate a wide 
range of jet flow conditions. 


One of the most popular techniques for computing supersonic jet flows 
is a space marching scheme based on solving the steady, parabolized Navier- 
Stokes (PNS) equations. The principle advantage of the PNS approach is its 
greater computational efficiency compared to that of methods which solve the 
full unsteady Navier-Stokes equations. The efficiency results from the fact 
that a solution can be obtained by a spatial marching approach in which the 
solution is advanced downstream from some specified initial condition. Thus, 
for steady supersonic flows, only a single marching sweep is needed. 

At the present time, a few such codes are available for predicting 
such flows, but they are limited in calculation capabilities. Dash and 
Wolf developed a two-dimensional (SCIPVIS (Ref. 10)) and a three-dimensional 
(SCIP3D Ref. 11)) parabolized Navier-Stokes code for analyzing propulsive 
jet mixing problems. SCIPVIS solves the mean flow equations for steady- 
state, two-dimensional compressible flows. This code can quantitively 
predict many of the details of the shock-cell structure, the turbulent 
mixing with an external stream, and the subsequent decay of the shock-cell 
strength as the result of shock/mixing-layer interactions. The SCIPVIS 
code can also give accurate predictions for underexpanded and overexpanded 
cases. SCIP3D solves the mean flow equations for steady-state, three- 
dimensional compressible flows. Wolf et al . (Ref. 12) have used SCIP3D 
to simulate an axisymetric supersonic jet problem. They found that, SCIP3D 
code overpredicted the shock-cell decay and spaciny as compared with the 
SCIPVIS result. These comparisons indicate that further work and modifi- 
cations in SCIP3D code are required to duplicate SCIPVIS result. They 
recommended the use of time-iterative procedure in plane to plane basis 
to replace their noniterative methodology which is proven to be very 
complex to deal with. They also suggested that the governing equation 
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should be formulated in generalized coordinates form, which will simplify the 
solution procedure for arbitrary exit shapes. 

In this study, an underexpanded supersonic jet will be used as a test 
case for developing the present three-dimensional code. The shock-cell 
structure in underexpanded, supersonic jets has been the subject of several 
experimental and theoretical studies in recent years (Ref. 13-17). An 
understanding of these structures is especially important to the field of 
aeroacoustics where shocks contribute significantly to jet noise (Ref. 13). A 
schematic of a typical flow field for an underexpanded, supersonic jet is 
given in fig. 1. The jet is characterized by repetitive shock cells whose 
strength and size are modified through turbulent mixing with the surrounding 
medi urn. 

The objective of the present study is the application and comparisons of 
two recently developed codes to simulate the shock-cell structure of the 
underexpanded supersonic jet problem. These codes were developed by modifying 
the 3D code developed by Newsome et al . (Ref. 1). Their code was used to 
solve the thin layer Navier-Stokes equations for a laminar, hypersonic 
afterbody flow problem. On the other hand, the present codes will be used in 
solving the Thin-layer Navier-Stokes (TLNS) and the Parabolized Navier-Stokes 
(PNS) equations for laminar or turbulent jet flow problems. Different options 
were added to the original code. These options are as follows: 

1. Either half (one symmetry; HP) or quarter (two symmetry plane; QP) of 
the cross-flow plane can be simulated. The QP option can be useful 
in calculating square or circular jets,. 

2. Van Leer's Flux Vector Splitting method or Rao's Flux Differencing 
technique can be used to calculate the flux in the subsonic and 
supersonic regions of the flow. 
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3. A limiter can be used in the shock region, thus eliminating 
overshoots or undershoots. 

4. The state variables at cell interface is constructed from either the 
primative or conservative variables. 

The present codes use an implicit, approximately-factored , upwind, flux 
splitting finite volume algorithm. The flux-splitting and upwind spatial 
differencing method for the convection terms, used in the present study, has 
several advantages over central difference schemes. This method has natural 
numerical dissipation and better stability property. Both codes are third 
order accuracy in the cross-plane directions and second order in the 
streamwise direction. These codes solve the equations in conservation form 
for a generalized coordinate system. The paper will also discuss some of the 
differences in computational requirements between the PNS and NS codes. 

GOVERNING EQUATIONS 

The TLNS equations can be written in the following generalized coordinate and 
conservative form: 


3Q 3E. 9^ 9G 

3t 35 3n 34 


0 


( 1 ) 


where 
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E = Vi + Vi + Vi * E i 


4 



F - R X (E 1 - E v ) ♦ R y (F, - F y ) ♦ R 2 ( Gj - G v ) = F, + 


G = T x E i + T y F i + T z G i ■ G i 


S,S,S,R,R,R»T,T andT are the components of the surface 
x y z x y* z x y z 

normals. These terms and the cell volume, Vol , are evaluated using the 
procedures described by Chakravarthy and Szema (Ref. 19). 

The inviscid (convection) flux vectors of the TLNS equations are given by 
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where the total energy is 
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The viscous (diffusion) flux vector are of the form 
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The Parabolized Navi er-Stokes (PNS) equations are obtained from the TLNS 
equations when the unsteady terms are omitted. If a space-marching procedure 
is used to solve these equations, the following assumptions should be enforced 
(Ref. 20): 

1. The streamwise velocity component is, everywhere, greater than zero 

2. The pressure gradient term in the streamwise direction (-|£) is either 
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omitted or treated with other technique to avoid a complex 
eigenval ue. 

In this study, the technique of Vingeron et. al . (Ref. 20) is adopted to 
supress the departure solutions associated with the elliptic behaviour of the 
equations. Vingeron et al . (Ref. 20) show that PNS equations is hyperbolic- 
parabolic provided that 


Ej = [pU, pull + S x wp, pvU + Sytop, puill + S z ojp, (Et + p)U] 


( 4 ) 
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where 


o> = 1 


oyM. 


1 + (y - 1) M 
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M_ > 1 

t. 


M < 1 


a is a safety factor for not taking into account the nonlinearity of the 
governing equations. 

COMPUTATIONAL METHOD 

In this section, the implicit upwind/ relaxation algorithm is used to solve the 
unsteady Navier-Stokes equation. The laws of conservation of mass, momentum, 
and energy over a volume, Vol , bounded by a surface S, can be expressed in 
integral form as 


It J vol « dV + A <* n)dS=0 (6) 

The tensor R can be written in terms of the Cartesian fluxes by 


A A A A A A 

B = (Ej - E v )i ♦ (F, - F v )j + (Gj - G v )k 

Associating the subscripts j , k , 1 with 5, n, C directions, a numerical 
approximation to Eq. (5) may be written in the following form: 
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where E, 6 and F are numerical fluxes at the boundary sides of the cell 

A 

and Q is the cell -average value for the numerical approximation to Q. 

A A A 

The flux vectors (E, F and G) are split according to the scheme of Van 
Leer (Ref. 21). These fluxed are split according to their contravariant Mach 
number 
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where 
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As an example, for supersonic flow in the £ direction 
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and for subsonic flow -1 < < 1 
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The flux at (n +1) is linearized as 


E n + 1 = £ n + a \ Q 


Now, the flux difference is written as 
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The conserved variables Q + , Q" obtained by an upwind biased one parameter 
fami ly 
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<(> = 0 


first order fully upwind 



second order fully upwind 


third order biased upwind 


However, to ensure Montone interpolation for the third order interpolation in 
the vicinity of a shock, a min mod limiter is used as follows: 


VQ = min mod (VQ,bAQ) 


AQ = min mod (AQ , bVQ) 


where b is a compression parameter, b = 



(13) 


It should be mentioned that the upwind split-flux difference procedure 

A A 

are only used for the inviscid convection parts of the flux vectors (F & G). 

A second order, central difference is used to represent the diffussion terms 
of SNS and PNS equations. 

In the present study, first order upwind method is used to evaluate the 
implicit (Jacobian part of equation (11)). This will result in great 
reduction of the computational steps and will not affect the steady-state 
accuracy which is controlled by the R.H.S. of the equation. Equation (6) can 
be written as 
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{ P V ^ , W^- > P^:}^ (P» “ 0,V,W,p}^ ^ 

These relations are equivalent to set the gradient of p, v, w and E and the 
value of u to zero along Sq line. 

3) Conditions at ? 

' max 

Factitious grid points are located along the side of s max line as shown 
in figure (3). For the quarter plane (QP) simulation, the values of the 
factitious points are 

T T 

lPf» U.r, Vp, Wp Pf) = ^P> o, - w » PS „ 

T T T T T n Sax’ n 

4) Condition at the centerline n Q or n = 0. (x = 0. and y = 0.) 

Factitious points are added across the n Q line as shown in figure (3). 
Their values are evaluated as follows: 


QP simulation 
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HP simulation 


(p f , u f , v f , w f , p f } T = (p - u, v, w, p}^ 
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5) Conditions on the Far-Field Boundary 


a) n 
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max 


The treatment of this boundary is based on Rieman invarients for a one- 
dimensional flow. 


Riemann invarients 
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Outflow Boundary, 0 >, 0 

K 

{P,u,v,w,p}^ = , u,v,w,p)3 + (V - V ) { 0 ,r ,r ,r ,o} T 

max a^ max 7 


This boundary is treated as the same as in 5-a where the surface normal 
r is replaced by s in the above calculations. These calculations are, 
only, used for the subsonic portion of the flow, whereas first order 
extrapolation was used for the supersonic part. 

RESULTS AND DISCUSSIONS 

In this section, the present three-dimensional codes are used to predict 
the shock-cell decay and flow characteristics of an underexpanded supersonic 
jet (1. < M < 2.) issued into still air. The three-dimensional results are, 
then, compared with the predicted results using a two-dimensional PNS code 
( SC I P V IS (Ref. 10) as well as the experimental data of Norum and Seiner (Ref. 
16) exhausting into nominally still air. It should be mentioned that PNS 
techniques assume that the velocity in the streamwise direction always has a 
non-negative value. Norum and Shearin (Ref. 22) investigate the shock 
structure and noise of supersonic jets in simulated flight to Mach 0.4. They 
show that the changes in shock strength and amplitude of broadband shock noise 
in flight to Mach 0.4 are insignificant. Therefore, all calculations 
presented in this paper are for an external stream at a Mach number of 0.05. 
Comparisons are presented for static pressure distributions measured along the 
jet centerline. 
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Two test cases were selected for comparison: 

1. Underexpanded, cold-air, supersonic jet -M = 2.0, p/p = 1.45. 

2. Underexpanded, cold-air, sonic jet - M = 1.0, p/p = 1.62 
In the present PNS calculations, 

1. The jet starts with a top-hat profile at x = 0. 

2. The local time-like iteration is activated at j-plane for n-iteration 

until the R.H.S. of the equation (15) is reduced to four-order of 
magnitude from its original value. 

3. The solution at j is used as an initial condition for the j + 1-plane 

which, in turn, reduce the number of iterations required to reach 

converged solution at J + 1-plane. 

4. Step 2 and 3 are repeated for all j-planes. Both test cases are 
converged within a number of interations less than 70. The PNS 
solutions are used as an initial condition for the TLNS solution as 
suggested by Walter et al . This approach will reduce the number of 
globel iteration required for the TLNS to reach a steady state 

sol ution. 

To get a good insight of the capability of each code using less computer 
resource (memory and time), we have selected a coarse grid in the T-direction 
to predict the shock-cell structure of the Mach 2 jet. The Quarter-Pi ane 
option will be sufficient in solving axisymmetric problems. A typical grid 
distribution is shown in figure 2 with 160x61x11 (JxKxL) grid. It is believed 
that the highest gradient (pressure, temperature, and streamwise velocity), in 
the i-direction occurs close to the sonic-line (sonic-surface for 3-D case). 

In this investigation, the location of the sonic-line is always less than 1.5 
jet-radii from the jet centerline. Based on these observations, a fine grid 
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was used for 1.5 jet-radii and stretched grid up to jet-radii. On the other 
hand, a uniform grid was used in the 5-direction. 

Case 1 - M = 2.0, P/P = 1.45 

In figure (4), the measured centerline static pressure illustrates the 
decaying shock structure that occurs as a result of the interaction of the 
shocks with the growing mixing layer. The jet is operated at a pressure ratio 
of 1.45 which corresponds to a fully expanded Mach number of 2.24 and is 
issued from a convergent-divergent nozzle with a design Mach number of 2. 
Figure (5) shows the comparison between the results predicted using PNS and 
TLNS code for the streamwise pressure variation along the jet centerline. The 
PNS predictions were completed in 240 CPU seconds compared to 8400 CPU seconds 
used by the TLNS code. This comparison shows that the PNS code is quite 
efficient, as compared to TLNS code, in predicting the shock-cell structure of 
this case. Figure (6) shows the comparison between the pressure and density 
contour, respectively, using PNS and TLNS code. In this figure, the PNS is 
able to give essentially the same result presented by the TLNS code. However, 
the TLNS code gives a better description for the compression and expansion 
(reflected and intercepting) shocks for the whole flow domain than the one 
given by the PNS solution which gives good description for the first 5 shock 
cel 1 s . 

Figure (7) shows the calculating centerline pressure using 160, 320 and 
640 grid points in the streamwise direction respectively. This comparison 
shows that increasing the grid refinement in the steamwise direction does not 

have a significant effect in the predictions of shock-cell spacing (number of 
shock-cell) however, it does show improvement in the strength and shape of the 

first three shock-cell (in particular the first reflected shock). 
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Finally, figure (8) shows the comparisons between 3-D and 2-D PNS 
solution with the measured streamwise pressure variation along the jet 
centerline. The measure static pressure distributions illustrate the decaying 
shock structure that occurs due to the interaction of the shocks with the 
growing mixing layer. Both the PNS solutions predict the pressure variations 
up through the first 4 shock cells and underestimate the rest of shock-cell 
decay. This was expected as the mixing layer and turbulent dissipation of the 
shock strength are not modeled in both codes. This problem can be resolved by 
using an adequate turbulence model. 

Case 2 - M = 1.0, P/P = 1.61 

The experimental results shown in figure (9) are for the underexpanded 
condition ( P j / Pa = 1.61) for a sonic nozzle. This pressure ratio corresponds 
to a fully expanded Mach number of 1.37, and this figure shows that the 
measured data appear to decay abruptly after the fourth shock cell. Although 
there are no flow-field data available to explain this uncharacteristic rapid 
decay, acoustic data obtained under similar conditions suggested that such 
phenomena may be caused by acoustic resonance. Figure (10) shows the results 
predicted with the 3-D and 2-D PNS codes compared with measured centerline 
pressures. Both calculations show reasonably good agreement up through the 
first 4 shock cells, although there is some disagreement in amplitude after 
the first shock cell. However, both calculations greatly underpredict the 
shock-cell decay after the fourth cell. It is believed that the flow did not 
reach steady state, and the measurement after the fourth shock cell is the 
representations of time average flows. It is believed that the disagreements 
with experiments are related to the limitations of the PNS methodology to 
handle the elliptic effects. 
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(See the Appendix for the fl ux-Jacobiaris J* , J* J* , and V„ 

1 u O L 

The implicit upwing/ relaxation algorithm of Newsom et. al (Ref. 1) is 
used to solve TLNS equations. This can be achieved through a series of 
alternative sweeps in the streamwise direction. For a forward 
sweep, AQ , is known 

U “ 1 5 K J I 

and AQ n ^ | , is set to zero. For a backward sweep, a. , , is known 

J + 1,K,I J + 1 ,K , I 

and AQ. 1 . , is set to zero. 

J “ i K 3 I 

Finally, equations (14) is approximately factored and can be written in 
the following compact form: 

A A 

(M + \ fjj) M" 1 (M + 6 c ||) AQ n + 1 = R.H..S. (15) 

where , 


M 


IVol 

At 


+ B 
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In the present code, the PNS equations are solved by retaining the time 
dependent terms and the steady state solution can be obtained by local 
iteration. This will allow the use of the same upwind flux-split procedure 
described in this section. Otherwise, very complicated flux-split procedure 
will be required because of the complex eigen-values and eigen-vectors of the 
steady state equations. At each cross-plane, local integration is used until 
the R.H.S. of equation (15) is reduced by 4-order of magnitude. Then, this 

solution is used as an initial condition for the next cross-sectional plane. 
This procedure is repeated for each cross-sectional plane cepta 5 max 
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THE BOUNDARY CONDITIONS 


Figure (2) shows a typical physical domain, used in the present study. 

The figure also shows the different boundary associated with external flow 
proplem. 

1) Inlet Flow Conditions at $ = 0. (z = 0) 

The static pressure, P. stream wise velocity, w, and temperature, are 
extracted form the measurement of Norum and Seiner (Ref. 14). The values of 
X- velocity, u, and y-velocity, v, are set to zero. Then, density p satisfies 
the equation of state 

p = P/R g e (16) 

and 

E t = ttt + h (lj2 + v2 + * 2) (17) 

which is used to calculate the total energy, E^. 

2) Conditions at 4 Q or 4 = 0 . (x = 0.) 

Factitious grid points are located along the side of 4 line as shown in 
figure (3). The points are denoted by subscript f, and their values are 
obtained from the following relations. 
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Figure (11) shows the comparisons between 3-D PNS and TLNS (1000 and 2000 
time steps) solutions for the sonic jet with pressure ratio of 1.61. The 
first 3 shock cells did not show any significant differences for the PNS and 
TLNS solutions. It is clear that the flow is unstable and will not reach a 
steady state. Figure (12) shows the density and pressure contours for the 
same results shown in figure (11). In this figure, the PNS and TLNS predict 
similar shock cell structure for the first 3 shock cells. For the TLNS 
solution, the first three shock cell seem to be quite stable and their 
structure is completely time independent. Close to the third shock cell a 
nearly normal shock is formed and it interacts with the basic structure. The 
flow beyond this point is completely unstable and vortex rings are formed on 
the jet boundary of the entire flow field. The same obsrervations were made 
by Matsuda et. al . (Ref. 17). Matsuda et. al . 17 investigated, 

experimentally and theoretically, the decay of the shock cell structure for a 
sonic jet with pressure ratio of 2. They took two kinds of Schieren 
photographs, a long time exposure ( 1/ 30s ) and and short time exposrue 
(1/lOOOs). With a short time exposure, they observed that the flow is 
turbulent with vortices (unstable). 

On the other hand, the long time exposure gives a rather regular and smooth 
shock cell structure. 

SUMMARY 

The development of PAB3D PNS and TLNS codes and their application to 
nonaxi symmetric jet mixing problem are described in this report. These codes 
were compared with the experimental data for underexpanded supersonic jets 
exhausting into a subsonic external stream.. The PNS code successfully 
predicts the shock cell structure of Mach 2.0 jet as compared with the TLNS 
code results. The PNS calculations were completed in less than 3 percent of 


21 



the TLNS computational time. However, PNS code fails to predict the flow 
field of strongly resonant jet (sonic case) beyond the first three shock 
cells. This was expected as the flow beyond this point was completely 

unstable and the only way to simulate this kind of flow is through the use of 

time dependent codes. For this kind of flow problem, the static or total 
pressure measurements should be accompanied by Schlieren photographs of a 
short and long time exposure or any method which can show the flow time 
dependency. 

In general, the PNS code can simulate efficiently, mildly underexpanded 
and fully supersonic flow problems, and can be used to set the initial 
conditions for a time dependent code. On the other hand, the TLNS can handle 

relatively complicated flow problems. The developments of these codes should 

be followed by a detailed validation study using nonaxi symmetric jet mixing 
data base. Right now, the TLNS code is currently being modified to include 
a multi-zone option which will give it the capability to simulate nozzle/ 
afterbody/ jet exhaust problems including vertical and horizontal tails. 
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APPENDIX 


This appendix decribes in detail the Jacobians of the PNS and TLNS governing 
equations for the inviscid and viscous terms 


1) Inviscid FI ux-Jacobians 


+ + 
3Et + 3Fr 3G: 

J 1 = TQ 5 J 2 = and = 


3Q 


3Q 


a) Supersonic M j> 


> 1 


J~ = 0 


3E 


J = 


1 

3Q 


< 1 


= 0 


j; = 


3Ei 

3Q 
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The above matrix was used for PNS and NS calculations. For PNS calculations. 


a) is evaluated using eq. (4) and J_ vanishes. On the other hand, 
set to value of 1 for the SNS calculations. 

where 

<j> = y 1 (u 2 + u 2 + w 2 ) , 0 = S u + S u + S w 
c. x y ^ 

D i = ( E t + p)/p 

b) Subsonic for SNS Solutions |M^| < 1 

(1.1) (1 - M 2 ) (a + P D 3 ) 

^ (1,2) = -D 4 D 2 u + S x D 5 

(1.3) = -D 4 0 2 v + S y D 5 

(1.4) = -D 4 D 2 w + S z D 5 
J* (1,5) = D 4 D 2 

J* (2,1) = J* (1,1) D 12 ±2s x p D 6 D 3 + t(s x D 7 - u) D 6 
0* (2,2) = 0 + (1.2)D 12 ^S x pD g D 2 u + y (1 - s 2 )D g 
^ (2,3) - J* (1,3)D 12 + 2s x pD 6 D 2 v - s x s y D 6 


(jo is 
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J* (2,4) = J* (1,4) D 12 +2s x pD 6 D 2 w - S x S z D g 
J* (2,5) = J± (1,5)D 12 ± 2s x pD 6 D 2 

(3,1) = J* (1,1) D 13 ± 2 s yP D 6 D 3 + Y(s y D ? - u)D 6 
J* (3,2) = J* (1,2)D 13 ?2s y pD g D 2 u - s x S y D g 
J* (3,3) = J* (1,3)D 13 + 2s y pD g D 2 v + Y D 6 (1 - s 2 ) 

J* (3,4) = J* (1,4) D 13 +2s yP D 6 D 2 w - s y s 2 D g 

J* (3,5) = J* (1,5)D 13 ± 2s yP D g D 2 

J* (4,1) = J* (1,1) ± 2s z pD g D 3 + y(s z D 2 - w) D g 

J* (4,2) = J* (1,2) D 15 +2s z pD 6 D 2 u - s z S x D g 

J* (4,3) = J* (1,3) D u +2 s z pD 6 D 2 v - s z s y D g 

J* (4,4) = J* (1,4) D 14 +2s z pD g D 2 w + y D 6 

J* (4,5) = J* (1,5) D 14 ± 2s z pD g D 2 

* 0 2<j> 

J (5,1) = J (1,1) D + (D D + D 

1 1 15 8 3 9 s (y ~1)p °11 
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J* (5,2) = J* (1,2) D 15 - (D 8 D 2u + s x O q + ^) D 


11 


(5.3) - (1,3) - (DgD^v + S yDg + p 

(5.4) = J* (1,4) D 15 - (D 8 D 2 w + s z D g ) D n 


(5,5) = J* (1,5) D 15 - D p D 2 D n 


° 2 - JLll ^r i ■ °3 - ° 2 <t^t - 4> 


d 5 = i < i ± ■ 


E* (1) 


’ "6 ip • D 7 = y 


Vttt <-5 ±a ^ 

Y - 1 


^ + * i 

q = £ “(1) D = ^ ~ D = — - 

U 11 L i U '* U 12 D„ 5 13 D, 


11 


11 


F ± 

n 1 (4) , n L i (5) 

D, „ = — - ---- and D„ v ' 


E. 


14 D 


11 


15 D 
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For the n direction, flux - Jacobian (J*) 


(o=l and 
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For the 5 direction FI ux-Jacobi an (j*), 

a>=l M -*■ M , 0 -*■ (9 , S T and Er ->• G* 

^ 11 s + t 

2) Viscous FI ux-Jacobians 

Neglecting cross-derivative term, the viscous flux vector in the n is given 
by: 


0 


F 


- m 
v ~ Vol 


M 

x n 


R 0 
y n 


R 0 
z n 


2 

+ R u 
2 

+ rS 

n 

2 

+ R w 


n 


R 2 4 > + U +4 R R (uv) + R R (vw) + R R ( wu ) _ + -7 - (a 2 ) 

Y n n 3 x y ' ' n y z v 'n z x ' '1 (y - 1) P ' '1 


where, Vol is the cell volume. 


0 =4(Ru + R v + R w ) , U =4- R 2 (u 2 ) + R 2 (v 2 ) + R 2 (w 2 ) 

n 3 x n yn z n' n 6 x v 'n yn z v/ n 


and , A (u 2 ) + (* 2 ) + (w 2 ) 

n 2 n n n 
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The viscous flux vector in the 5 direction can be constructed by replacing the 
following variables in the above formulas, 


F 

v 


G 


v 


R + T n + c and R = [R x> 



T 
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The viscous fl ux-Jacobian is given as follows: 
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Figure 1. Typical flow field for underexpanded supersonic jet. 
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Figure 2. Typical grid distribution for the physical domain. 






PRESSURE (ATM) 




Figure 5. Comparison of 3-D PNS and 3-D NS predicted centerline pressures for 
underexpanded supersonic jet. 
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A) PNS SOLUTION 



DENSITY CONTOUR 


CONTOUR INTERVAL - .11900E+00 


B) NS SOLUTION 



DENSITY CONTOUR 


CONTOUR INTERVAL = .11942E+0Q 


Figure 6a. Comparison of 3-D and 3-D PNS predicted density contours for 
underexpanded supersonic jet. 
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A) PNS SOLUTION 



PRESSURE CONTOUR 


CONTOUR INTERVAL - .37381 E+04 


B) NS SOLUTION 



PRESSURE CONTOUR 


CONTOUR INTERVAL - .37929E+04 


Figure 6b. Comparison of 3-D and 3-D PNS predicted pressure contours for 
underexpanded supersonic jet. 
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Figure 7. Comparison of 3-D PNS predicted centerline pressures for 
underexpanded supersonic jet. 
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Figure 8. Comparison of predicted (2-D and 3-D PNS) and measured centerline 
pressure for underexpanded supersonic jet. 
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Figure 9. Measured plume centerline pressure for underexpanded sonic jet. 
(Data from ref. 14) 
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Figure 10. Comparison of predicted (2-D and 3-D PNS) and measured centerline 
pressure for underexpanded sonic jet,, 
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Figure 12a. Comparison of 3-D and 3-D PNS predicted density contours for 
underexpanded sonic jet. 
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Figure 12b. Comparison of 3-D and 3-D PNS predicted pressure contours for 
underexpanded sonic jet. 
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